Lista 2 - Econometria 2

Questão 1: Asset Pricing

(a)

## DADOS DE FAMA-FRENCH:

# 1. Importação dos Dados
ff_raw <- read_csv("25_Portfolios_5x5.csv", skip = 15)

# 2. Limpeza e Tratamento Mensal
ff_monthly <- ff_raw %>%
  rename(Date = `...1`) %>%
  filter(!is.na(Date), !str_detect(Date, "Copyright")) %>%
  mutate(Date = ymd(paste0(Date, "01"))) %>%
  mutate(across(-Date, as.numeric))

# Carregar e preparar os dados de CPI
cpi_raw <- read_csv("CPIAUCNS.csv")

# Cálculo do componente inflacionário
cpi_monthly <- cpi_raw %>%
  rename(CPI = CPIAUCNS) %>%
  mutate(Date = ymd(observation_date)) %>%
  arrange(Date) %>%
  mutate(Inflation_pct = (CPI / lag(CPI) - 1) * 100) %>%
  select(Date, Inflation_pct)

# Adicionando o CPI ao dataframe
ff_monthly_real <- ff_monthly %>%
  left_join(cpi_monthly, by = "Date") %>%
  filter(!is.na(Inflation_pct))

# 3. Conversão para Trimestral (compondo os retornos com o componente inflacionário)
ff_quarterly_real <- ff_monthly_real %>%
  mutate(YearQuarter = as.yearqtr(Date)) %>%
  group_by(YearQuarter) %>%
  summarise(
    across(
      .cols = -c(Date, Inflation_pct),
      .fns = ~ ((prod(1 + .x / 100, na.rm = TRUE) / prod(1 + Inflation_pct / 100, na.rm = TRUE)) - 1) * 100
    ),comp_inflation_quarterly = (prod(1 + Inflation_pct / 100, na.rm = TRUE) - 1) * 100
  ) %>%
  ungroup()

head(ff_quarterly_real)

## # A tibble: 6 × 27
##   YearQuarter `SMALL LoBM` `ME1 BM2` `ME1 BM3` `ME1 BM4` `SMALL HiBM` `ME2 BM1`
##   <yearqtr>          <dbl>     <dbl>     <dbl>     <dbl>        <dbl>     <dbl>
## 1 1947 Q1             40.8      26.9      51.6      71.1         114.      41.3
## 2 1947 Q2             66.6      46.8     107.      120.          205.      78.4
## 3 1947 Q3             54.7      71.6     119.      217.          374.     112. 
## 4 1947 Q4             32.8      92.5     148.      273.          429.     125. 
## 5 1948 Q1             88.4     104.      162.      295.          499.     138. 
## 6 1948 Q2             73.9     149.      178.      302.          511.     162. 
## # ℹ 20 more variables: `ME2 BM2` <dbl>, `ME2 BM3` <dbl>, `ME2 BM4` <dbl>,
## #   `ME2 BM5` <dbl>, `ME3 BM1` <dbl>, `ME3 BM2` <dbl>, `ME3 BM3` <dbl>,
## #   `ME3 BM4` <dbl>, `ME3 BM5` <dbl>, `ME4 BM1` <dbl>, `ME4 BM2` <dbl>,
## #   `ME4 BM3` <dbl>, `ME4 BM4` <dbl>, `ME4 BM5` <dbl>, `BIG LoBM` <dbl>,
## #   `ME5 BM2` <dbl>, `ME5 BM3` <dbl>, `ME5 BM4` <dbl>, `BIG HiBM` <dbl>,
## #   comp_inflation_quarterly <dbl>

## DADOS DE CONSUMO:

# 1. Importação do Consumo Nominal (Tabela 2.3.5)

# Carregar os cabeçalhos
header_nominal <- read_csv("Table 2.3.5. Personal Consumption Expenditures by Major Type of Product.csv",
  skip = 3,
  n_max = 2,
  col_names = FALSE
)

# Carregar os dados
data_nominal <- read_csv(
  "Table 2.3.5. Personal Consumption Expenditures by Major Type of Product.csv",
  skip = 5,
  col_names = FALSE,
  col_types = cols(.default = "c")
)

# Criar nomes de colunas corretos
year_row_nom <- header_nominal[1, -c(1, 2)]
quarter_row_nom <- header_nominal[2, -c(1, 2)]
new_column_names_nom <- paste0(zoo::na.locf(unlist(year_row_nom)), "_", unlist(quarter_row_nom))
colnames(data_nominal) <- c("Line", "Description", new_column_names_nom)

# Isolar a série de consumo nominal de não duráveis (Linha 8)
nondurables_nominal <- data_nominal %>%
  filter(Line == "8") %>%
  pivot_longer(
    cols = contains("_Q"),
    names_to = "Quarter_str",
    values_to = "Nondurables_Nominal"
  ) %>%
  mutate(
    YearQuarter = as.yearqtr(Quarter_str, format = "%Y_Q%q"),
    Nondurables_Nominal = as.numeric(Nondurables_Nominal)
  ) %>%
  select(YearQuarter, Nondurables_Nominal)


# 2. Importação do Índice de Preços (Tabela 2.3.4)

# Carregar os cabeçalhos
header_prices <- read_csv("Table 2.3.4. Price Indexes for Personal Consumption Expenditures by Major Type of Product.csv",
  skip = 3,
  n_max = 2,
  col_names = FALSE
)

# Carregar os dados
data_prices <- read_csv(
  "Table 2.3.4. Price Indexes for Personal Consumption Expenditures by Major Type of Product.csv",
  skip = 5,
  col_names = FALSE,
  col_types = cols(.default = "c")
)

# Criar nomes de colunas corretos
year_row_price <- header_prices[1, -c(1, 2)]
quarter_row_price <- header_prices[2, -c(1, 2)]
new_column_names_price <- paste0(zoo::na.locf(unlist(year_row_price)), "_", unlist(quarter_row_price))
colnames(data_prices) <- c("Line", "Description", new_column_names_price)

# Isolar a série do índice de preços de não duráveis (Linha 8)
price_index_nondurables <- data_prices %>%
  filter(Line == "8") %>%
  pivot_longer(
    cols = contains("_Q"),
    names_to = "Quarter_str",
    values_to = "Price_Index"
  ) %>%
  mutate(
    YearQuarter = as.yearqtr(Quarter_str, format = "%Y_Q%q"),
    Price_Index = as.numeric(Price_Index)
  ) %>%
  select(YearQuarter, Price_Index)


# 3. Calculando o Consumo Real

# Unir as duas tabelas e deflacionar os valores nominais
nondurables_quarterly <- nondurables_nominal %>%
  inner_join(price_index_nondurables, by = "YearQuarter") %>%
  mutate(
    # A fórmula para deflacionar é (Valor Nominal / Índice de Preços) * 100
    Nondurables_Real = (Nondurables_Nominal / Price_Index) * 100
  ) %>%
  select(YearQuarter, Nondurables_Real)

cat("--- Amostra do Consumo Real de Não-Duráveis Calculado ---\n")

## --- Amostra do Consumo Real de Não-Duráveis Calculado ---

head(nondurables_quarterly)

## # A tibble: 6 × 2
##   YearQuarter Nondurables_Real
##   <yearqtr>              <dbl>
## 1 1947 Q1                 502.
## 2 1947 Q2                 510.
## 3 1947 Q3                 512.
## 4 1947 Q4                 506.
## 5 1948 Q1                 507.
## 6 1948 Q2                 513.

# 4. Importação e Adequação do Dataframe de População (Tabela 2.1)
header_pop <- read_csv("Table 2.1. Personal Income and Its Disposition.csv",
  skip = 3,
  n_max = 2,
  col_names = FALSE
)

data_pop <- read_csv(
  "Table 2.1. Personal Income and Its Disposition.csv",
  skip = 5,
  col_names = FALSE,
  col_types = cols(.default = "c")
)

# Criando os nomes de coluna corretos
year_row_pop <- header_pop[1, -c(1, 2)]
quarter_row_pop <- header_pop[2, -c(1, 2)]
new_column_names_pop <- paste0(zoo::na.locf(unlist(year_row_pop)), "_", unlist(quarter_row_pop))
colnames(data_pop) <- c("Line", "Description", new_column_names_pop)

# 5. Limpeza e Transformação para a série de População
population_quarterly <- data_pop %>%
  filter(Line == 40) %>%
  pivot_longer(
    cols = contains("_Q"),
    names_to = "Quarter_str",
    values_to = "Population_Thousands"
  ) %>%
  mutate(
    YearQuarter = as.yearqtr(Quarter_str, format = "%Y_Q%q"),
    Population_Thousands = as.numeric(Population_Thousands)
  ) %>%
  select(YearQuarter, Population_Thousands)

# 6. Unindo os dataframes de consumo e população
consumption_full <- nondurables_quarterly %>%
  left_join(population_quarterly, by = "YearQuarter")

# Verificando se há NAs na nova coluna
cat("\nValores Faltantes (NA) na coluna de População:", sum(is.na(consumption_full$Population_Thousands)), "\n")

## 
## Valores Faltantes (NA) na coluna de População: 0

# 7. Criando a coluna de consumo per-capita
consumption_full <- consumption_full %>%
  mutate(Nondurables_Real_PC = (Nondurables_Real * 1000000000) / (Population_Thousands * 1000))

summary(consumption_full$Nondurables_Real_PC)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3446    4496    5850    6030    7686   10256

head(consumption_full)

## # A tibble: 6 × 4
##   YearQuarter Nondurables_Real Population_Thousands Nondurables_Real_PC
##   <yearqtr>              <dbl>                <dbl>               <dbl>
## 1 1947 Q1                 502.               143143               3505.
## 2 1947 Q2                 510.               143790               3548.
## 3 1947 Q3                 512.               144449               3546.
## 4 1947 Q4                 506.               145122               3487.
## 5 1948 Q1                 507.               145709               3481.
## 6 1948 Q2                 513.               146289               3504.

# 8. Calculando o Crescimento do Consumo Real Per-Capita
consumption_final <- consumption_full %>%
  arrange(YearQuarter) %>%
  mutate(
    consumption_growth_pc = log(Nondurables_Real_PC) - lag(log(Nondurables_Real_PC))
  )

head(consumption_final)

## # A tibble: 6 × 5
##   YearQuarter Nondurables_Real Population_Thousands Nondurables_Real_PC
##   <yearqtr>              <dbl>                <dbl>               <dbl>
## 1 1947 Q1                 502.               143143               3505.
## 2 1947 Q2                 510.               143790               3548.
## 3 1947 Q3                 512.               144449               3546.
## 4 1947 Q4                 506.               145122               3487.
## 5 1948 Q1                 507.               145709               3481.
## 6 1948 Q2                 513.               146289               3504.
## # ℹ 1 more variable: consumption_growth_pc <dbl>

## DADOS DOS FATORES FAMA-FRENCH E TAXA LIVRE DE RISCO (RF):

# 1. Importação dos Dados
factors_raw <- read_csv("F-F_Research_Data_Factors.csv", skip = 3)

# 2. Limpeza e Tratamento Mensal
factors_monthly <- factors_raw %>%
  rename(Date = `...1`) %>%
  filter(!is.na(Date)) %>%
  mutate(Date = ymd(paste0(Date, "01"))) %>%
  mutate(across(-Date, as.numeric))

# 3. Conversão da RF para Trimestral (compondo a taxa)
rf_quarterly <- factors_monthly %>%
  mutate(YearQuarter = as.yearqtr(Date)) %>%
  group_by(YearQuarter) %>%
  summarise(
    RF_quarterly = (prod(1 + RF / 100, na.rm = TRUE) - 1) * 100
  ) %>%
  ungroup()


## JUNÇÃO FINAL E CÁLCULO DO SDF

# 1. Unir as bases de dados (Fama-French, Consumo e Taxa Livre de Risco)
final_data <- ff_quarterly_real %>%
  left_join(consumption_final, by = "YearQuarter") %>%
  left_join(rf_quarterly, by = "YearQuarter")

# 2. Criar o Dataframe Final de Análise
analysis_df <- final_data %>%
  # 2.1 Calcular o excesso de retorno (prêmio de risco) para cada portfólio
  mutate(across(
    .cols = starts_with("SMALL") | starts_with("ME") | starts_with("BIG"),
    .fns = ~ .x - RF_quarterly,
    .names = "{.col}_excess"
  )) %>%
  # 2.2 Calcular o SDF
  mutate(
    gross_consumption_growth = exp(consumption_growth_pc),
    SDF = 0.99 * (gross_consumption_growth ^ -2)
  ) %>%
  # 2.3 Remover a primeira linha que tem NA no SDF (devido ao lag do crescimento)
  filter(!is.na(SDF))

# 3. Calcular o Prêmio de Risco Estimado pelo Modelo
mean_sdf <- mean(analysis_df$SDF, na.rm = TRUE)

model_premiums <- analysis_df %>%
  summarise(across(
    .cols = ends_with("_excess"),
    .fns = ~ -cov(SDF, .x, use = "complete.obs") / mean_sdf,
    .names = "model_{.col}"
  ))

model_premiums_long <- model_premiums %>%
  pivot_longer(cols = everything(), 
               names_to = "Portfolio", 
               values_to = "Model_Estimated_Premium")

cat("--- Prêmios de Risco Trimestrais Estimados pelo Modelo\n")

## --- Prêmios de Risco Trimestrais Estimados pelo Modelo

print(model_premiums_long, n = 25)

## # A tibble: 25 × 2
##    Portfolio               Model_Estimated_Premium
##    <chr>                                     <dbl>
##  1 model_SMALL LoBM_excess                  490.  
##  2 model_SMALL HiBM_excess                  497.  
##  3 model_ME1 BM2_excess                       2.63
##  4 model_ME1 BM3_excess                      63.1 
##  5 model_ME1 BM4_excess                     182.  
##  6 model_ME2 BM1_excess                    1177.  
##  7 model_ME2 BM2_excess                    2260.  
##  8 model_ME2 BM3_excess                     631.  
##  9 model_ME2 BM4_excess                    2110.  
## 10 model_ME2 BM5_excess                     421.  
## 11 model_ME3 BM1_excess                   11420.  
## 12 model_ME3 BM2_excess                   23340.  
## 13 model_ME3 BM3_excess                   14701.  
## 14 model_ME3 BM4_excess                   16046.  
## 15 model_ME3 BM5_excess                    3249.  
## 16 model_ME4 BM1_excess                  414332.  
## 17 model_ME4 BM2_excess                  222608.  
## 18 model_ME4 BM3_excess                  145208.  
## 19 model_ME4 BM4_excess                  122116.  
## 20 model_ME4 BM5_excess                   66316.  
## 21 model_ME5 BM2_excess               154703655.  
## 22 model_ME5 BM3_excess                17136945.  
## 23 model_ME5 BM4_excess                30761523.  
## 24 model_BIG LoBM_excess             1210029838.  
## 25 model_BIG HiBM_excess               -8191002.

(b)

# 1. Tratando os Dados dos Fatores (análogo ao tratamento para os 25 portifólios da letra (a))
factors_monthly_raw <- read_csv("F-F_Research_Data_Factors.csv", skip = 3)

factors_monthly <- factors_monthly_raw %>%
  rename(Date = `...1`) %>%
  filter(!is.na(Date), !str_detect(Date, "Copyright")) %>%
  mutate(Date = ymd(paste0(Date, "01"))) %>%
  mutate(across(-Date, as.numeric))

factors_monthly_real <- factors_monthly %>%
  left_join(cpi_monthly, by = "Date") %>%
  filter(!is.na(Inflation_pct))

factors_quarterly_real <- factors_monthly_real %>%
  mutate(YearQuarter = as.yearqtr(Date)) %>%
  group_by(YearQuarter) %>%
  summarise(
    across(
      .cols = c(`Mkt-RF`, RF),
      .fns = ~ ((prod(1 + .x / 100, na.rm = TRUE) / prod(1 + Inflation_pct / 100, na.rm = TRUE)) - 1) * 100,
      .names = "{.col}_real"
    ),
    `Mkt-RF` = (prod(1 + `Mkt-RF` / 100, na.rm = TRUE) - 1) * 100,
    SMB = (prod(1 + SMB / 100, na.rm = TRUE) - 1) * 100,
    HML = (prod(1 + HML / 100, na.rm = TRUE) - 1) * 100
  ) %>%
  ungroup()

cat("--- Fatores Fama-French Trimestrais (Reais e Nominais) ---\n")

## --- Fatores Fama-French Trimestrais (Reais e Nominais) ---

head(factors_quarterly_real)

## # A tibble: 6 × 6
##   YearQuarter `Mkt-RF_real` RF_real `Mkt-RF`    SMB   HML
##   <yearqtr>           <dbl>   <dbl>    <dbl>  <dbl> <dbl>
## 1 1947 Q1            -4.53   -1.77    -2.75  -0.911 0.751
## 2 1947 Q2            -1.17   -0.365   -0.717 -7.42  0.525
## 3 1947 Q3            -2.66   -4.23     1.77   3.45  4.30 
## 4 1947 Q4             1.69   -1.51     3.45  -3.62  4.97 
## 5 1948 Q1            -0.688   0.230   -0.688  0.898 6.02 
## 6 1948 Q2             7.89   -2.66    11.1   -2.55  5.72

# 2. Unindo os Retornos dos Portfólios com o Fator de Mercado
capm_data <- analysis_df %>%
  select(YearQuarter, ends_with("_excess")) %>%
  inner_join(factors_quarterly_real %>% select(YearQuarter, `Mkt-RF`), by = "YearQuarter")

# 3. Estimando o CAPM para cada um dos 25 Portfólios
capm_results <- capm_data %>%
  pivot_longer(
    cols = ends_with("_excess"),
    names_to = "Portfolio",
    values_to = "Excess_Return"
  ) %>%
  # Agrupa por portfólio para rodar uma regressão para cada um
  group_by(Portfolio) %>%
  # tidyr::nest() cria um "sub-dataframe" para cada portfólio
  nest() %>%
  # purrr::map() aplica a função de regressão (lm) a cada sub-dataframe
  mutate(
    model = map(data, ~ lm(Excess_Return ~ `Mkt-RF`, data = .x)),
    # broom::tidy() extrai os coeficientes (alpha e beta) de cada modelo de forma organizada
    tidied = map(model, tidy)
  ) %>%
  # Expande a lista de coeficientes de volta para o formato de tabela
  unnest(tidied)

cat("\n\n--- Coeficientes do CAPM (Alpha e Beta) para cada Portfólio ---\n")

## 
## 
## --- Coeficientes do CAPM (Alpha e Beta) para cada Portfólio ---

print(capm_results)

## # A tibble: 50 × 8
## # Groups:   Portfolio [25]
##    Portfolio         data     model  term  estimate std.error statistic  p.value
##    <chr>             <list>   <list> <chr>    <dbl>     <dbl>     <dbl>    <dbl>
##  1 SMALL LoBM_excess <tibble> <lm>   (Int…  127810.    12113.    10.6   1.88e-22
##  2 SMALL LoBM_excess <tibble> <lm>   `Mkt…     217.     1273.     0.170 8.65e- 1
##  3 SMALL HiBM_excess <tibble> <lm>   (Int…  375921.    28416.    13.2   5.31e-32
##  4 SMALL HiBM_excess <tibble> <lm>   `Mkt…    8194.     2986.     2.74  6.41e- 3
##  5 ME1 BM2_excess    <tibble> <lm>   (Int…   69114.     6924.     9.98  1.54e-20
##  6 ME1 BM2_excess    <tibble> <lm>   `Mkt…    1460.      727.     2.01  4.56e- 2
##  7 ME1 BM3_excess    <tibble> <lm>   (Int…   79808.     7416.    10.8   3.61e-23
##  8 ME1 BM3_excess    <tibble> <lm>   `Mkt…    1647.      779.     2.11  3.54e- 2
##  9 ME1 BM4_excess    <tibble> <lm>   (Int…  141324.    11955.    11.8   7.02e-27
## 10 ME1 BM4_excess    <tibble> <lm>   `Mkt…    2961.     1256.     2.36  1.90e- 2
## # ℹ 40 more rows

# 4. Calculando o Prêmio de Risco Previsto pelo CAPM
# Calculamos a média do prêmio de risco de mercado histórico
avg_market_premium <- mean(capm_data$`Mkt-RF`, na.rm = TRUE)

# Filtramos para pegar apenas os Betas (o coeficiente de `Mkt-RF`) de cada portfólio
betas <- capm_results %>%
  filter(term == "`Mkt-RF`") %>%
  select(Portfolio, beta = estimate)

# Calculamos o prêmio de risco de cada portfólio: Beta * Média do Prêmio de Mercado
capm_risk_premiums <- betas %>%
  mutate(CAPM_Premium = beta * avg_market_premium)

cat("\n\n--- Prêmio de Risco Trimestral Estimado pelo CAPM\n")

## 
## 
## --- Prêmio de Risco Trimestral Estimado pelo CAPM

print(capm_risk_premiums, n = 25)

## # A tibble: 25 × 3
## # Groups:   Portfolio [25]
##    Portfolio                 beta CAPM_Premium
##    <chr>                    <dbl>        <dbl>
##  1 SMALL LoBM_excess         217.         487.
##  2 SMALL HiBM_excess        8194.       18390.
##  3 ME1 BM2_excess           1460.        3278.
##  4 ME1 BM3_excess           1647.        3696.
##  5 ME1 BM4_excess           2961.        6645.
##  6 ME2 BM1_excess          35473.       79616.
##  7 ME2 BM2_excess          46909.      105283.
##  8 ME2 BM3_excess          27898.       62613.
##  9 ME2 BM4_excess          30516.       68490.
## 10 ME2 BM5_excess          18057.       40526.
## 11 ME3 BM1_excess         477641.     1072009.
## 12 ME3 BM2_excess         598170.     1342521.
## 13 ME3 BM3_excess         236064.      529817.
## 14 ME3 BM4_excess         221389.      496881.
## 15 ME3 BM5_excess          77299.      173488.
## 16 ME4 BM1_excess       14184253.    31834846.
## 17 ME4 BM2_excess        4805607.    10785607.
## 18 ME4 BM3_excess        2732268.     6132247.
## 19 ME4 BM4_excess        2320861.     5208892.
## 20 ME4 BM5_excess        1379408.     3095915.
## 21 ME5 BM2_excess     2616557520.  5872548148.
## 22 ME5 BM3_excess      359205227.   806192860.
## 23 ME5 BM4_excess      986420947.  2213902985.
## 24 BIG LoBM_excess   17557041403. 39404664404.
## 25 BIG HiBM_excess     214126499.   480581132.

(c)

# 1. Preparando os Dados para a Regressão de 3 Fatores
# Unir os retornos reais dos portfólios com os três fatores Fama-French.
three_factor_data <- analysis_df %>%
  select(YearQuarter, ends_with("_excess")) %>%
  inner_join(
    factors_quarterly_real %>% select(YearQuarter, `Mkt-RF_real`, SMB, HML), 
    by = "YearQuarter"
  )

cat("--- Amostra dos Dados para o Modelo de 3 Fatores ---\n")

## --- Amostra dos Dados para o Modelo de 3 Fatores ---

head(three_factor_data)

## # A tibble: 6 × 29
##   YearQuarter `SMALL LoBM_excess` `SMALL HiBM_excess` `ME1 BM2_excess`
##   <yearqtr>                 <dbl>               <dbl>            <dbl>
## 1 1947 Q2                   66.5                 205.             46.7
## 2 1947 Q3                   54.6                 374.             71.5
## 3 1947 Q4                   32.6                 428.             92.3
## 4 1948 Q1                   88.2                 499.            104. 
## 5 1948 Q2                   73.7                 511.            148. 
## 6 1948 Q3                    2.34                320.             50.5
## # ℹ 25 more variables: `ME1 BM3_excess` <dbl>, `ME1 BM4_excess` <dbl>,
## #   `ME2 BM1_excess` <dbl>, `ME2 BM2_excess` <dbl>, `ME2 BM3_excess` <dbl>,
## #   `ME2 BM4_excess` <dbl>, `ME2 BM5_excess` <dbl>, `ME3 BM1_excess` <dbl>,
## #   `ME3 BM2_excess` <dbl>, `ME3 BM3_excess` <dbl>, `ME3 BM4_excess` <dbl>,
## #   `ME3 BM5_excess` <dbl>, `ME4 BM1_excess` <dbl>, `ME4 BM2_excess` <dbl>,
## #   `ME4 BM3_excess` <dbl>, `ME4 BM4_excess` <dbl>, `ME4 BM5_excess` <dbl>,
## #   `ME5 BM2_excess` <dbl>, `ME5 BM3_excess` <dbl>, `ME5 BM4_excess` <dbl>, …

# 2. Estimando o Modelo de 3 Fatores para cada Portfólio (análogo à letra (b))
three_factor_results <- three_factor_data %>%
  pivot_longer(
    cols = ends_with("_excess"),
    names_to = "Portfolio",
    values_to = "Excess_Return"
  ) %>%
  group_by(Portfolio) %>%
  nest() %>%
  mutate(
    model = map(data, ~ lm(Excess_Return ~ `Mkt-RF_real` + SMB + HML, data = .x)),
    tidied = map(model, tidy)
  ) %>%
  unnest(tidied)

cat("\n\n--- Coeficientes do Modelo de 3 Fatores para cada Portfólio ---\n")

## 
## 
## --- Coeficientes do Modelo de 3 Fatores para cada Portfólio ---

print(three_factor_results, n = 12)

## # A tibble: 100 × 8
## # Groups:   Portfolio [25]
##    Portfolio         data     model  term  estimate std.error statistic  p.value
##    <chr>             <list>   <list> <chr>    <dbl>     <dbl>     <dbl>    <dbl>
##  1 SMALL LoBM_excess <tibble> <lm>   (Int…  133129.    12060.    11.0   4.17e-24
##  2 SMALL LoBM_excess <tibble> <lm>   `Mkt…     154.     1315.     0.117 9.07e- 1
##  3 SMALL LoBM_excess <tibble> <lm>   SMB     -1082.     2173.    -0.498 6.19e- 1
##  4 SMALL LoBM_excess <tibble> <lm>   HML     -4456.     2013.    -2.21  2.76e- 2
##  5 SMALL HiBM_excess <tibble> <lm>   (Int…  378767.    28449.    13.3   2.82e-32
##  6 SMALL HiBM_excess <tibble> <lm>   `Mkt…    8756.     3103.     2.82  5.08e- 3
##  7 SMALL HiBM_excess <tibble> <lm>   SMB      1112.     5127.     0.217 8.28e- 1
##  8 SMALL HiBM_excess <tibble> <lm>   HML      2867.     4749.     0.604 5.47e- 1
##  9 ME1 BM2_excess    <tibble> <lm>   (Int…   70548.     6943.    10.2   4.09e-21
## 10 ME1 BM2_excess    <tibble> <lm>   `Mkt…    1404.      757.     1.85  6.47e- 2
## 11 ME1 BM2_excess    <tibble> <lm>   SMB       500.     1251.     0.400 6.90e- 1
## 12 ME1 BM2_excess    <tibble> <lm>   HML      -318.     1159.    -0.274 7.84e- 1
## # ℹ 88 more rows

# 3. Calculando o Prêmio de Risco Previsto pelo Modelo de 3 Fatores
# Calculamos a média histórica de cada um dos três fatores
avg_factor_premiums <- three_factor_data %>%
  summarise(
    avg_mkt_real = mean(`Mkt-RF_real`, na.rm = TRUE),
    avg_smb = mean(SMB, na.rm = TRUE),
    avg_hml = mean(HML, na.rm = TRUE)
  )

# Extraímos os betas de cada portfólio para cada fator
betas_3f <- three_factor_results %>%
  filter(term %in% c("`Mkt-RF_real`", "SMB", "HML")) %>%
  select(Portfolio, term, estimate) %>%
  pivot_wider(names_from = term, values_from = estimate) %>%
  rename(
    beta_mkt = `\`Mkt-RF_real\``,
    beta_smb = SMB,
    beta_hml = HML
  )

# Calculamos o prêmio de risco para cada portfólio
three_factor_premiums <- betas_3f %>%
  mutate(
    Three_Factor_Premium = 
      beta_mkt * avg_factor_premiums$avg_mkt_real +
      beta_smb * avg_factor_premiums$avg_smb +
      beta_hml * avg_factor_premiums$avg_hml
  )

cat("\n\n--- Prêmio de Risco Trimestral Estimado pelo Modelo de 3 Fatores\n")

## 
## 
## --- Prêmio de Risco Trimestral Estimado pelo Modelo de 3 Fatores

print(three_factor_premiums, n = 25)

## # A tibble: 25 × 5
## # Groups:   Portfolio [25]
##    Portfolio             beta_mkt      beta_smb    beta_hml Three_Factor_Premium
##    <chr>                    <dbl>         <dbl>       <dbl>                <dbl>
##  1 SMALL LoBM_excess         154.        -1082.     -4.46e3               -4832.
##  2 SMALL HiBM_excess        8756.         1112.      2.87e3               15544.
##  3 ME1 BM2_excess           1404.          500.     -3.18e2                1844.
##  4 ME1 BM3_excess           1619.          386.     -7.22e2                1677.
##  5 ME1 BM4_excess           2073.         6802.      4.13e2                6388.
##  6 ME2 BM1_excess          36715.        -4253.     -7.11e3               41571.
##  7 ME2 BM2_excess          45972.         4310.      2.59e3               68188.
##  8 ME2 BM3_excess          30694.       -12397.      2.63e3               39490.
##  9 ME2 BM4_excess          31687.        -3243.      7.59e3               50092.
## 10 ME2 BM5_excess          20864.       -11356.      1.29e4               36873.
## 11 ME3 BM1_excess         484325.       -31923.      6.30e4              719604.
## 12 ME3 BM2_excess         587048.        65336.      1.50e5              995153.
## 13 ME3 BM3_excess         243748.       -51657.      3.23e4              346590.
## 14 ME3 BM4_excess         224658.        -2379.      3.49e4              345264.
## 15 ME3 BM5_excess          94616.       -85674.      5.66e4              149613.
## 16 ME4 BM1_excess       13808157.      1505782.      1.59e6            21410852.
## 17 ME4 BM2_excess        5002452.     -1267525.      6.04e5             6957196.
## 18 ME4 BM3_excess        2850086.      -340073.      1.60e6             5416120.
## 19 ME4 BM4_excess        2565810.     -1424794.      6.56e5             3568838.
## 20 ME4 BM5_excess        1632878.     -1239227.      1.56e6             3285212.
## 21 ME5 BM2_excess     2665886771.   -387086610.      8.43e8          4371596994.
## 22 ME5 BM3_excess      456205123.   -541643157.     -1.75e8           205289163.
## 23 ME5 BM4_excess     1104354255.   -633080922.      3.57e8          1603921830.
## 24 BIG LoBM_excess   20678555373. -20940905374.     -4.96e9         13990816776.
## 25 BIG HiBM_excess     316906469.   -546665652.      1.39e8           330948133.

2.

ipca<-ipeadata("PRECOS12_IPCA12")%>%filter(date>=as.Date("2000-01-01"),date<=as.Date("2025-01-01"))%>%arrange(date)
saveRDS(ipca,"ipca_precos12_2000_2025.rds");write.csv(ipca,"ipca_precos12_2000_2025.csv",row.names=FALSE)

a)

ipca<-readRDS("ipca_precos12_2000_2025.rds");s<-ts(ipca$value,start=c(2000,1),frequency=12)
summary(ipca$value);c(mean=mean(ipca$value),sd=sd(ipca$value),min=min(ipca$value),max=max(ipca$value))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1598    2574    3498    3858    5093    7112

##     mean       sd      min      max 
## 3858.085 1577.834 1598.410 7111.860

ggplotly(autoplot(s)+theme_classic()+ggtitle("IPCA - Índice (2000-2025)"))%>%layout(plot_bgcolor="white",paper_bgcolor="white",yaxis=list(rangemode="tozero"))

ggplotly(autoplot(log(s))+theme_classic()+ggtitle("Log do IPCA"))%>%layout(plot_bgcolor="white",paper_bgcolor="white",yaxis=list(rangemode="tozero"))

ggplotly(autoplot(diff(log(s)))+theme_classic()+ggtitle("Δlog(IPCA)"))%>%layout(plot_bgcolor="white",paper_bgcolor="white",yaxis=list(rangemode="tozero"))

Pode-se observar uma notável tendência de crescimento do IPCA tanto em nível quanto em log, enquanto o comportamento dos dados em diferença do log apresentam um comportamento mais oscilatório, sem tendência clara de crescimento e se aproximando um pouco mais de um ruído branco.

adf.test(s);pp.test(s)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  s
## Dickey-Fuller = -0.43059, Lag order = 6, p-value = 0.9846
## alternative hypothesis: stationary

## 
##  Phillips-Perron Unit Root Test
## 
## data:  s
## Dickey-Fuller Z(alpha) = -0.29147, Truncation lag parameter = 5,
## p-value = 0.99
## alternative hypothesis: stationary

Pode-se observar que os testes não rejeitam a hipótese nula de raiz unitária, dado o elevado p-valor dos testes, o que indica uma iminente tendência estocástica e não-estacionariedade, conforme indicado pelos gráficos.

adf.test(log(s));pp.test(log(s))

## 
##  Augmented Dickey-Fuller Test
## 
## data:  log(s)
## Dickey-Fuller = -2.3204, Lag order = 6, p-value = 0.4413
## alternative hypothesis: stationary

## 
##  Phillips-Perron Unit Root Test
## 
## data:  log(s)
## Dickey-Fuller Z(alpha) = -7.4022, Truncation lag parameter = 5, p-value
## = 0.6953
## alternative hypothesis: stationary

Assim como em nível, em log se observa um p-valor elevado, com os testes não rejeitando a hipótese nula de raiz unitária, indicando novamente tendência estocástica e não-estacionariedade

adf.test(diff(log(s)));pp.test(diff(log(s)))

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(log(s))
## Dickey-Fuller = -5.813, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff(log(s))
## Dickey-Fuller Z(alpha) = -118, Truncation lag parameter = 5, p-value =
## 0.01
## alternative hypothesis: stationary

Para os dados em diferença do log, observa-se um p-valor consideravelmente menor, rejeitando a hipótese nula de raiz unitária e indicando estacionariedade.

d<-diff(log(s));a<-acf(d,plot=FALSE,na.action=na.pass);p<-pacf(d,plot=FALSE,na.action=na.pass);ci<-qnorm(0.975)/sqrt(a$n.used)
xa<-as.numeric(a$lag)*12;ya<-as.numeric(a$acf);xp<-as.numeric(p$lag)*12;yp<-as.numeric(p$acf)
plot_ly(x=xa,y=ya,type="bar")%>%layout(title="ACF Δlog(IPCA)",plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Defasagem (meses)"),yaxis=list(title="Correlação",rangemode="tozero"),shapes=list(list(type="line",xref="x",yref="y",x0=min(xa),x1=max(xa),y0=ci,y1=ci),list(type="line",xref="x",yref="y",x0=min(xa),x1=max(xa),y0=-ci,y1=-ci)))

plot_ly(x=xp,y=yp,type="bar")%>%layout(title="PACF Δlog(IPCA)",plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Defasagem (meses)"),yaxis=list(title="Correlação",rangemode="tozero"),shapes=list(list(type="line",xref="x",yref="y",x0=min(xp),x1=max(xp),y0=ci,y1=ci),list(type="line",xref="x",yref="y",x0=min(xp),x1=max(xp),y0=-ci,y1=-ci)))

Enquanto na PACF não se observa nenhum lag significativo, com todos os lags após o lag nulo se mantendo dentro do intervalo de confiança, pode-se assumir que as partes de média móvel sejam 0 tanto no aspecto não-sazonal quanto no sazonal, com possivelmente \(q=Q=0\).

Contudo, analisando a ACF, é possível observar um decaimento gradual após o lag 0 e além disso um comportamento de crescimento após um ano, indicando um possível fator AR tanto sazonal quanto não-sazonal, indicando possivelmente que \(p=P=1\).

Dessa forma, o modelo que tende a ser o mais adequado para o IPCA é o \(SARIMA(1,1,0)(1,0,0)_{12}\).

b)

library(forecast);library(plotly)
y<-log(s)
fit1<-Arima(y,order=c(1,1,0),seasonal=list(order=c(1,0,0),period=12),include.drift=FALSE,method="ML")
fit2<-Arima(y,order=c(1,1,0),seasonal=list(order=c(0,0,1),period=12),include.drift=FALSE,method="ML")
summary(fit1);summary(fit2)

## Series: y 
## ARIMA(1,1,0)(1,0,0)[12] 
## 
## Coefficients:
##          ar1    sar1
##       0.8060  0.2269
## s.e.  0.0375  0.0656
## 
## sigma^2 = 1.06e-05:  log likelihood = 1295.33
## AIC=-2584.65   AICc=-2584.57   BIC=-2573.54
## 
## Training set error measures:
##                        ME        RMSE         MAE         MPE       MAPE
## Training set 0.0007735025 0.003238879 0.002341851 0.009619096 0.02868748
##                    MASE       ACF1
## Training set 0.03888154 -0.1301465

## Series: y 
## ARIMA(1,1,0)(0,0,1)[12] 
## 
## Coefficients:
##          ar1    sma1
##       0.8232  0.1775
## s.e.  0.0336  0.0559
## 
## sigma^2 = 1.069e-05:  log likelihood = 1294.07
## AIC=-2582.14   AICc=-2582.05   BIC=-2571.02
## 
## Training set error measures:
##                        ME        RMSE         MAE         MPE       MAPE
## Training set 0.0007725472 0.003253154 0.002351729 0.009602986 0.02882324
##                    MASE      ACF1
## Training set 0.03904554 -0.138668

Estimados os dois modelos, podemos observar valores muito próximos nos testes de AIC e BIC, mas o modelo com termo AR sazonal apresenta um valor substancialmente menor em ambos os testes, sendo esse o modelo mais indicado com base nos dois testes.

plres<-function(f,t){
  r<-na.omit(residuals(f));ix<-seq_along(r)
  a<-acf(r,plot=FALSE);p<-pacf(r,plot=FALSE);ci<-qnorm(0.975)/sqrt(a$n.used)
  xa<-as.numeric(a$lag)*12;ya<-as.numeric(a$acf);xp<-as.numeric(p$lag)*12;yp<-as.numeric(p$acf)
  rt<-plot_ly(x=ix,y=as.numeric(r),type="scatter",mode="lines")%>%layout(title=paste0("Resíduos - ",t),plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Índice temporal"),yaxis=list(title="Resíduo",rangemode="tozero"))
  ac<-plot_ly(x=xa,y=ya,type="bar")%>%layout(title=paste0("ACF dos resíduos - ",t),plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Defasagem (meses)"),yaxis=list(title="Correlação",rangemode="tozero"),shapes=list(list(type="line",xref="x",yref="y",x0=min(xa),x1=max(xa),y0=ci,y1=ci),list(type="line",xref="x",yref="y",x0=min(xa),x1=max(xa),y0=-ci,y1=-ci)))
  pc<-plot_ly(x=xp,y=yp,type="bar")%>%layout(title=paste0("PACF dos resíduos - ",t),plot_bgcolor="white",paper_bgcolor="white",xaxis=list(title="Defasagem (meses)"),yaxis=list(title="Correlação",rangemode="tozero"),shapes=list(list(type="line",xref="x",yref="y",x0=min(xp),x1=max(xp),y0=ci,y1=ci),list(type="line",xref="x",yref="y",x0=min(xp),x1=max(xp),y0=-ci,y1=-ci)))
  list(rt=rt,ac=ac,pc=pc)
}
g1<-plres(fit1,"SARIMA(1,1,0)(1,0,0)[12]")
g2<-plres(fit2,"SARIMA(1,1,0)(0,0,1)[12]")
g1$rt;g2$rt

g1$ac;g2$ac

g1$pc;g2$pc

Dado que os dois modelos tem coeficientes próximos, é esperado um desempenho similar dos gráficos de resíduos, algo que de fato ocorre, com ambos aproximando-se de um ruído branco. Em ambos os modelos se observa um decaimento após o lag 0 com os lags \(1\) e \(2\) sendo negativos e fora do intervalo de confiança por uma diferença mínima além de um lag \(15\) significativo, o que poderia talvez indicar a necessidade de estimação de mais parâmetros para um ajuste mais fino, mas que perderia em parcimônia com a estimação de novos parâmetros. O mesmo vale para os valores dos lags \(1\) e \(2\) fora do intervalo de confiança da PACF em ambos os modelos.

lb<-function(f,k)c(LB24=Box.test(residuals(f),type="Ljung-Box",lag=24,fitdf=k)$p.value,
                   LB36=Box.test(residuals(f),type="Ljung-Box",lag=36,fitdf=k)$p.value)
rbind("(1,1,0)(1,0,0)[12]"=lb(fit1,2),
      "(1,1,0)(0,0,1)[12]"=lb(fit2,2))

##                          LB24      LB36
## (1,1,0)(1,0,0)[12] 0.07489526 0.2695822
## (1,1,0)(0,0,1)[12] 0.06336111 0.2285484

Os testes de Ljung-box para os lags \(24\) e \(36\) apresentaram p-valores superiores a \(5\%\), o que indica ausênsia de correlação residual, e então que os modelos estão capturando adequadamente a estrutura temporal da série, sem restar um padrão sistemático claro dos erros aparentemente.